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Abstract 

An analytical soluble model based on a Continuous Time Random Walk (CTRW) scheme for the 
adsorption-desorption processes at interfaces, called bulk-mediated surface diffusion, is presented. 
The time evolution of the effective probability distribution width on the surface is calculated 
and analyzed within an anomalous diffusion framework. The asymptotic behavior for large times 
shows a sub-diffusive regime for the effective surface diffusion but, depending on the observed range 
of time, other regimes may be obtained. Montecarlo simulations show excellent agreement with 
analytical results. As an important byproduct of the indicated approach, we present the evaluation 
of the time for the first visit to the surface. 

PACS numbers: 
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I. INTRODUCTION 

The dynamics of adsorbed molecules is a fundamental issue in interface science 111 141 and is 

nn 

crucial to a number of emerging technologies LjJ • Its role is central to phenomena as diverse 
as foam relaxation (f| and the evolution of blood protein deposit Q|. Recently, the mechanism 
called bulk-mediated surface diffusion has been identified and explored. Its importance in 
relaxing homogeneous surface density perturbations is experimentally well established. This 
mechanism arises at interfaces separating a liquid bulk phase and a second phase which 
may be either solid, liquid, or gaseous. Whenever the adsorbed species is soluble in the 
liquid bulk, adsorption-desorption processes occur continuously. These processes generate a 
surface displacement because molecules desorb, undergo Fickian diffusion in the bulk liquid, 
and are then re-adsorbed elsewhere. When this process is repeated many times, it results 
in an effective diffusion of a molecule on the surface. Bichuk and OShaughnessy jsj have 
claimed that this effective surface diffusion has anomalous super-diffusive characteristics 
when certain range of time is considered. 

Dynamical processes that display anomalous diffusion have been 

characterized by a non linear time dependence of the mean square displacement of the walker 
for long times; that is (r 2 (t)) ~ t e with t / l (remember that e = 1 corresponds to normal 
diffusion) since (r 2 ) is the usual estimator of the square width of the probability distribution 
at time t. Hence, for anomalous diffusion we have that the probability distribution width 
grows faster (slower) for e > 1 (e < 1) than it does for normal diffusion. 

In this paper we present an analytical soluble model for the adsorption-desorption pro- 
cesses based on a Continuous Time Random Walk (CTRW) scheme. We calculate the 
evolution with time of the square width of the effective probability distribution on the sur- 
face and show that, for a given range of time, this square width growth as t e where e depends 
on the values of the adsorption and diffusion parameters. 



II. THE ADSORPTION-DESORPTION MODEL 



Let us start with the problem of a particle making a random walk in the semi- infinite 
cubic lattice. The position of the walker is defined by the vector r whose components are 
denoted by the integer numbers n,m,l > 1 corresponding to the directions x, y and z 
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respectively. The displacements in the x and y directions are unbounded. In the z direction 
the particle can move from I — 1 to infinity. 

The probability that the walker is in (n, m, Z) at time t, given that it was at (0, 0, Zo) at to, 
P{n, m, Z; t\n, m, Zo; to) = P(n, m, I; t), satisfies the following set of coupled master equations 

P{n,m,l;t) = jP(n, m, 2; t) — SP(n, m, 1; t), 1 = 1 

P(n, m, 2; t) = a[P(n - 1, m, 2; *) + P(n + 1, m, 2; t) - 2P(n, m, 2; t)] + 
/3[P(n, m - 1, 2; t) + P(n, m + 1, 2; t) - 2P(n, m, 2; *)] + 
7P(n, m, 3; t) + 6P(n, m, 1; t) - 27P(n, m, 2; t), Z = 2 

P(n, m, /; t) = a[P(n - 1, m, Z; t) + P(n + 1, m, /; i) - 2P(n, m, I; t)] + 
/3[P(n, m-l,l;t) + P(n, m+l,l;t)- 2P(n, m, /; *)] + 
7[P(n,m,/ + l;t) + P(n,m,/-l;t)-2P(n,m,/;t)], Z>3 (1) 

where a, f3 and 7 are the transition probabilities per unit time in the x,y and z directions 
respectively, and 5 is the desorption probability per unit time from the boundary plane 
defined by z = 1. 

Taking the Fourier transform with respect to the x and y variables and the Laplace 
transform with respect to the time t in the above equations, we obtain 

sG(k x , k y , 1; s) - P(k x , k y , 1, t = 0) = 7G(A; X , 2; s) - 5G(k x , k y , 1; s), Z = 1 
sG[k x , ky, 2; s) P(k x , ky, 2, t — 0) = A(k x , ky)G(k x , ky, 2; s) + 8G\k x , ky, 1; s) + 

^G(k x , k y , 3; s) - 2^G(k x , k y , 2;s), 1 = 2 
sG{k x , ky, Z; s) — P(k x , ky, l,t — 0) = A(A;.j., ky)G{k x , ky, Z; s) + 

7[G(fca., A; y , Z - l;s) + 

^(^^Z + ljs) -2G(A; :c ,A; 2/ ,Z;s)], Z>3 (2) 

We have used the following definitions 

e" s ' ^ e k * n+kym P(n,m,l;t\n,m,l ;to)dt 

n,m=— 00 

00 

= L[ ^ e fc ^ m P(n,m,Z;t|n,m,Z ;to)], (3) 

n,m=— 00 

where L indicates the Laplace transform of the quantity within the brackets, and 

fc„) = 2a[Cos(A; :c ) - 1] + 2/3[Cos(A; 2/ ) - 1]. (4) 
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It is possible to write Eq. (J2J) in matrix form as 

[uI-H]G = 5 Uo = I Uo , (5) 
where the square matrix G has components 

Gu = [G[k x , k y , I; s\n, m, l ; £„]], (6) 
In Eq. (|5ji . I is the identity matrix and H is a three-diagonal matrix with the following form 

( -5 7 ...\ 

5 C 7 ••• 

fl"= 7 C 7 ••• 

7 C 7 ••■ 

V / 

and C is defined as 

C = -2 1 + A(k x ,k y ). (7) 

In order to find the solution to the Eq. (jSJ), we decompose the H matrix in the following 

way 

H = A(k x , k y )i + H + Hi + H 2 , (8) 

where 

/-27 7 
7 -2 7 7 
i?o = 7 —27 7 
7 -2 7 7 

V / 



{Hxh = Ai 



1 ifi=j = 1 
otherwise 



with 



(#2)^ = A 2 



1 if i = 1 and j = 2 
otherwise 



Ai = -5 - [-27 + A(k x , ky)}, 
A 2 = 5-7. 



(9) 



Defining 

G° = [(s-{A{k x ,k y ))I-H )}-\ 

G 1 = [(s - {A{k x , k y ))I -H - Hi)]' 1 , (10) 
observing that a formal solution of Eq. (JHJ) is 

G — [si — H]" 1 , (11) 
and by reiterating the Dyson formula, we can show that 

Gu„ = Gi„ + P^£, (12) 



The form for G® 1q can be obtained by conventional methods j3| as 

= L(e 2 ^(%_ i0 |(2 7i ) _ / b . + . o| (2 7 t))) fi , (14) 

where I n {2^ft) is the first modified Bessel function of order n. The above expression points 
out that the Laplace transform is evaluated at the argument u = [s — A(k x , k y )]. 

Once the general expression for G is obtained, we can find the probability that a particle 
is on the plane at site (m,n,l = 1) at time t given it was at (0,0,/ = 1) at t = 0. This 
probability may be obtained using the inverse Laplace transform in s and the inverse Fourier 
transform on k x , k y of the matrix element G\\. 



A direct measurable experimental magnitude 
bution at time t over the plane z = 1 



is the variance of the probability distri- 



< r 2 (t) > p i ane , (15) 

which measures the spreading of particles over this plane. Once P(m,n,l = l;t|0,0,/o = 
1; t — 0) is known, the variance is calculated as: 

oo 

<r 2 (t)> plane = p ( m ,n,l = l;t|0,0,Z = l;t )(m 2 + n 2 ). (16) 



m,n=— oo 



Here, we have used symmetry properties for the diffusion along the x and y axes, that is 
< x {t) >=< y{t) >= 0. 



The Laplace transform of the variance < r 2 (s) > p i an e= L(< r 2 (t) > p i a ne) can be found as 
follows 

d 2 d 2 ~ 

< r 2 (s) > p i ane = -[-^ + Qjp][Gn] \k x =k y =o (17) 

x y 

By using Eqs. (JTTJ) to (fTTjl . < r 2 {u) > p i an e turns out to be the ratio of two complicated 
functions of s 

9/ x N(s) , x 

< r 2 (s) > p l ane = (18) 

where 



N(s) = [45 7 (« + /?)][s 2 ( v / s(4 7 + s) - 3s) + 3 7 2 ( v / s(4 7 + s) - 3s) - 2 7 3 ], (19) 

and 



D(s) = y/s(4rf + s)( 7 s(2 7 + s - v / s(4 7 + s)) +<5( 7 ( v / s(4 7 + s) -3s) + s( v / s(4 7 + s)-s))) 2 . 

(20) 

It is important to remark that the conservation of particles in the plane is not satisfied. 

If we denote with P p i an e{t) the probability that the particle is in the plane z — 1 at time 
t, it can be shown that the Laplace transform of the magnitude is 

6 / \ ( 7 (2 7 + s- y/s(4 7 + s)) 

( 7 s(2 7 + s - v/s(4 7 + s)) + 5( 7 (v/s(4 7 + s) - 3s) + s(v/s(4 7 + s) - s))) 

(21) 

In order to test the theoretical results for the < r 2 {t) > p i ane and P p i an e(t), we have 
performed Montecarlo simulations for the adsorption-desorption processes obtaining an ex- 
cellent agreement in both cases. Figures [I] and 121 show the variance and the probability that 
the particle is on the plane Z — 1 clS cl function of £, for three different values of 5. 

The asymptotic behavior for large t of both, < r 2 (t) > p i an e and P p i an e(t), can be obtained 
by means of Tauberian theorems |l2| as 



Pplane{t) ~> ^ ^[1/2]^ ^ ^ ^ 23 ' ) 

From Eq. (|22j) we recognize an asymptotic sub-diffusive regime which is shown in Fig. El 

When we choose a different range of time in order to fit < r 2 (t) > p i an e as t e we find that 
e depends on the values of 5 for fixed a, (3 and 7. Figure |U shows this dependence for a wide 
range of values of the desorption rate S. 




FIG. 1: < r 2 (t) > vs time, for three different values of 5. Circles correspond to simulations 
with 5 = 0.01, while the solid line is the related theoretical curve. Squares are for simulations with 
5 = 0.1, and the dashed line is the theoretical curve. Triangles are for simulations with 5 = 1.0, and 
the dotted line is the corresponding theoretical curve. In all figures we have used a = (3 = 7 = 1.0, 
and the number of realizations is 10 6 . 

III. THE CTRW SCHEME 

In this section we sumarize the most important results of the lattice CTRW approach. 
Let \l/(r, r';t) the waiting time density governing a single transition being denned as the 
probability density that a transition from r' to f occurs at time between t and t + dt given 
the walker arrives at r at t — 0. The transition time between different sites is assumed much 
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FIG. 2: Probability P(z = l;i) vs time, for three different values of 5. Circles correspond to 
simulations with 5 = 0.01, while the solid line is the theoretical curve. Squares are for simulations 
with 5 = 0.1, and the dashed line is the theoretical curve. Triangles are for simulations with 
8 = 1.0, and the dotted line is the corresponding theoretical curve. 

shorter than the time spent at each lattice point; these "instantaneous transitions" allow us 
to prove probability conservation as indicated below. 

We assumed translational invariance, that is ^(r, r';t) = \l/(f — r'\ t) for all r and r', 
hence the conditional probability to find the particle at site r at time t given it was at r = 
at t — 0, P(r, t\f— 0, t — 0) can be found in Laplace- Fourier variables k and s as 

P(k,s\r = 0;t = 0) = -U— 1- ^ (g) (24) 

l-w(M) s 
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FIG. 3: < r 2 (t) > vs time, for three different values of 5, and compared with a fit using y = at € . 
Circles correspond to theory with 5 = 0.01, and the solid line to the fit with a = 0.63 ± 0.03, 
e = 0.663 ± 0.007. Squares are for theory with S = 0.1, while the solid line is the fit with 
a = 0.435 ± 0.007, e = 0.535 ± 0.005. The insert shows the related results for simulations. Again 
circles correspond to 5 = 0.01, and the solid line to the fit with a = 0.64±0.04 and e = 0.661±0.008; 
while squares are for 5 = 0.1, and the solid line is the fit with a = 0.47 ± 0.01, e = 0.531 ± 0.004. 

where ^(k, s) is the Fourier-Laplace transform of \l/(r, t) and ^(s) = ^{k = 0, s) is the 
Laplace transform of the total transition probability. 

The normalization of P(f,t\f — 0, t — 0) for all t is obtained from Eq. P4"|) in a direct 
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FIG. 4: e vs 5, for time in the range t = to t = 500. Note that for 5 — ► oo, e approaches the value 
0.5, while for small values of 5, e becomes ^> 1. 

way as 

^P(r]t\r = 6;t = 0) = L _1 (P(fc = 0, s|r = 0, t = 0)) = L _1 (-) = 1. (25) 

Normalization implies that the walker is somewhere in the (3D) lattice, but if we evaluate 
the probability that the walker is in a subspace of this lattice we will lose the normalization 
condition. 

We now consider, as done previously j*>], the possibility that the successive visits of the 
walker to the plane z = 1, as discussed in Section |Hj may be viewed as an CTRW over that 
plane. We must define a suitable waiting time density for single transitions, ^ p i ane (r,r';t) 
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where r and r' are two dimensional vectors on the plane, i.e. r— (m,n) and r' = (m',n'). 
A transition between r and r' must be done with no visit to the plane within the interval 
(0, £), in order to be consistent with a single transition. If we remember that this waiting 
time density is the probability density that the walker arrives at f between t and t + dt 
having arrived at r' at t — 0, it is obvious that this transition is not instantaneous because 
the "flying time" across the bulk cannot be neglected (see below). In order to be consistent 

— * 

with the CTRW theory we may consider that the particle remains at site r' during this 
time and then perform a jump to the site r. In this way, the probability over the plane is 
conserved and transitions become instantaneous. 

Now we build up the waiting time density for this "single transition" in the plane z — 1, 
taking into account the above remarks. We note that if the walker has arrived at (m',n f , 1) 
at t — 0, the probability density to desorb from the plane per unit time around t' by a jump 
to (to', n', 2) is ~ S exp(— St'). We define q[(m, n, 2), t; (to', n', 2), t'] as the probability of 
finding the walker at (to, n, 2) at time t given it was in (to', n', 2) at time t' without visiting 
the plane z — 1 in the interval (t',t). The probability density to reach (m,n, 1) for the first 
time between t and t + dt given the walker was in (to', n', 2) at time t' can be expressed as 

/[(to, n, l),t; (to', n', 2), t')] = >yq[(m, n, 2), t; (to', n', 2),t'}. (26) 

Finally the density to reach (to, n, 1) for first time, per unit of time around t given that 
the walker was in (m',n f , 1) at t — without visiting the plane in the interval (0,t) is: 

f[(m,n,l),t;(m',n',l),0]= I [ 7 g[(m, n, 2), t; (to', n', 2), t']]S exp(-5lf)dlf. (27) 

Jo 

The function /[(to, n, 1), t; (to', n', 1), t'] plays the role of the waiting time density 
^piane(r, r'\ t) where r = (m,n) and r' = (m',n'). Since we are assuming translational 
invariance in the x and y directions, the function q depends only on (to — to') and (n — n') . 
Selecting (0, 0) as the starting point we obtain 

^pianeim, n,t) = [ [ 79 (to, 7i, 2), t; (0, 0, 2), t'}5 exp(-5t')dt' . (28) 
Jo 

The function q((m, n, 2), t; (0, 0, 2), t') is obtained by means of the method of images assum- 
ing an absorbent plane in z — 1 

g(TO,n,2),t;(0,0,2),t') = [exp(-2at)J ro (2at)] 

[exp(-2/3t)I n (2/3t)] 

[exp(-27t)/ (2 7 t) - exp(-2 7 t/ 2 (2 7 t)], (29) 
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where Ij is the modified Bessel function of order j. Equations (|2%j) and (j29|) allow us to 
build the "normalized probability" of the CTRW in the plane z = 1 

P™(r, t\r = 0, t = 0) = L- 1 ^- 1 ( 1 _ ) f 1 ~ **™( s ) ) ( 30 , 



where 



^«ne(fc, S) = 7(^)(S + 2 7 + A(fc x , - yJ( S + 2 7 - A(^, fcj,)* - (2 7 )2))_, (31) 

and ^pioneCs) = ^pjanel^ = 0, s) . Here r = (to, n); = and the function v4(/c x , fcj,) is 

defined by Eq. (jlj). Eq. (|3~Tf shows a "coupled" waiting time density with a divergent time 
first moment, that is < t >—>■ oo. 

Expressions for the Fourier-Laplace transform of p norm (f, t\r — 0, i = 0) can be obtained 
from Eqs. (pjl and ([3Tj) . 

1(2 7 ( 5 + 5) - 5(5 + 2 7 - ^3(3 + 47))) 



s|r = 0,t = 0) 



2 7 (s + 5) - <S(s + 2 7 + A(k x , k y ) - ^J(s + 2 1 + A(k x ,k y )y-(2 1 y) 

(32) 



The corresponding variance of p norm in Laplace space is 

< r\s) >™= 25(a + ^±±±=^^^ ay /sJ^j), (33) 

2 7 s + o y s(s + 4 7 J — s 

and the asymptotic behavior for large t is 

< r 2 {t) (a + /?)£. (34) 

This result, i.e. normal diffusion, is due to the coupled character of the waiting time density 
Eq. (|31jl. and its infinite time first moment. 

As an important byproduct of the above approach, we present the evaluation of the 
probability of the first return to the plane z — 1. If a walker, initially at the point (0, 0, 1), 
desorbs and begins an excursion across the bulk, the probability density to return, for the 
first time, to the plane z = 1 between t and t + dt, f re t(t), is given by 

f ret (t) = L'^pian^s)] = L- 1 [^ one (^ = 0, s)]. (35) 

From the Eq. (|3~T|) we obtain the following expression the Laplace transform of the first 
return density 

6 w 1 



^planet*) = ( — )(— -)(2 7 + 8 ~ y/u( S + 4<y)). (36) 

2 7 s + o 
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FIG. 5: Distribution function for return times f re t(t) vs time for 5 = 0.5. The squares correspond 
to simulations while the solid line is for theory. 

We have made a numerical inverse transform of this result by using a numerical program 
and have compared this result with Montecarlo simulations finding an excellent agreement. 
See Fig. 

Finally, from Eq. ((HE)) it is possible to obtain two important results. Firstly, a walker is 
certain to return to the plane 

POO 

/ fret(t) dt = V p l ane (s = 0) = 1. (37) 
J 

Secondly, the asymptotic (long time) behavior of the first return density is 

/„*(*) - ( ^FI7 )4, (38) 
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indicating that the mean first time to return to the plane is infinite. 



IV. CONCLUSIONS 

We presented in this paper an analytical model for the adsorption-desorption processes 
from a boundary plane in a semi-infinite cubic lattice. We studied the effective diffusion 
of molecules on this plane interface and calculated the evolution of the square width of the 
probability distribution. This square width can be fitted with a power law of t whose expo- 
nent changes with the range of time considered and depends on the values of the adsorption 
and diffusion parameters. In this sense, effective anomalous super-diffusions reported by 
Bichuk and OShaughnessy (..§] may be understood. However, in the asymptotic long time 
regime this square width always behaves as t^. It is important to observe that the lack of 
probability conservation over the interface must be taken into account if a "genuine" CTRW 
on the plane is considered as was pointed out in Section ITTT1 As a byproduct of our approach, 
we obtained the evaluation of the probability of the first return to the planar interface. We 
performed Montecarlo simulations of the adsorption-desorption process obtaining excellent 
agreement with the model predictions. 

This work is part of a research project on bulk mediated diffusion on surfaces. Here 
we have discussed the case of infinite bulk while in we investigated the finite (in the 
direction normal to the surface) case where, among other aspects, we have found that an 
"optimal" number of layers exists, that produces the faster growth of < r 2 (t) >. In addition, 
we have also investigated the case of non-Markovian desorption process |19( where we have 
found an interesting oscillatory behavior. This research offers a more or less complete view 
of the theoretical description for the problem of bulk mediated diffusion on a surface. 
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